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Abstract 



We have studied the stability of the ferromagnetic state in the innnite-t/ 
Hubbard model on a square lattice by approximate diagonalization of finite 
lattices using the density matrix renormalization group technique. By study- 
ing lattices with up to 100 sites, we have found the ferromagnetic state to be 
stable below the hole density of 5 C = 0.22. Beyond 5 C , the total spin of the 

ground state decreased gradually to zero with increasing hole density. 
PACS: 71.27.+a,75.10.-b 
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The origin of many unusual electronic properties of high T c superconductors can be traced 
to strong electron-electron repulsion in the CuO planes. The Hubbard model is the simplest 
description of such repulsive interactions. Because of its simplicity the Hubbard model plays 
a role in many-body problems similar to that of the Ising model in phase transition problems. 
However, the Hubbard model is still very difficult to analyze. After forty years of research, 
we are still unsure of even its most basic features |]||. In this letter, we focus on the 
infinite-?/ limit to begin looking for unusual behaviors suggested recently There are 

several reasons for studying this limit. First of all, the antiferromagnetic state at half filling 
in the large U limit is incompatible with the motion of holes in the metallic phase |J. It is 
interesting to learn about the spin background preferred by the motion of the holes without 
the complication of the antiferromagnetic interaction. A well understood infinite-?/ limit 
also provides a starting point for systematic expansion in t/U. Furthermore, the ground 
state of the infinite-?/ Hubbard model at small doping is believed to be ferromagnetic. This 
provides a mechanism for itinerant ferromagnetism . But it is controversial whether there 
is a finite range of hole density where the ground state is ferromagnetic. This letter provides 
strong evidence that such a finite region indeed exists for the square lattice. We used the 
recently developed density matrix renormalization group (DMRG) method to compute the 
critical hole doping |7j]. 

The investigation of the infinite-?/ Hubbard model has a long history. The earliest rig- 
orous result is due to Nagaoka 0, and independently Thouless 0, showing that in the case 
of one hole, the ground state on a bipartite lattice is the ferromagnetic state (also known as 



2 



Nagaoka state), where all the spins are aligned in the same direction. Since then, progress 
on this difficult problem has been slow [p~0| p~6[] . Recently, it has been shown |I?,|IB| that for 
two holes the Nagaoka state is not the ground state. However, the proposed two-hole trial 



state [17fl has essentially local ferromagnetic correlation. Shastry et al. considered [19[ the 



instability of one spin flip of the Nagaoka state at finite hole density. They shown that the 
Nagaoka state is unstable when the hole density exceeds 5 C = 0.49. This result has been 
improved |20| to yield 5 C = 0.41. The single spin flip state has also been studied by von der 



Linden and Edwards ]21] using a more general trial wave function. They shown that the 
Nagaoka state is unstable against a single spin flip for 5 > 0.29. 

By comparing the high temperature expansion coefficients of the infinite-f/ Hubbard 



model with that of a free spinless fermion Hamiltonian, Yedidia [22] conjectured that the 
transition to the non-ferromagnetic state occurs at 5 ~ 3/11. The high temperature ex- 



pansion has been extended to higher order by Putikka et al. ] 23|| . When the free energy is 
extrapolated to zero temperature, their calculation suggests 8 C = 0. 



An extensive exact diagonalization investigation p4| , |25| has been carried out using Lanc- 
zos method, which is limited to small clusters. A very large finite-size effect has been ob- 
served. On a square lattice with periodic boundary conditions, the Nagaoka state is stabilized 
for the close-shell configurations when the number of holes is 1, 5, 9, .... At other hole fillings 
the Nagaoka state tends to be destabilized on small lattices because the energy change from 
one shell to the next is too large. Because of this, the ground-state magnetization oscillates 
with the number of holes 124" . 
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We have studied the stability of the Nagaoka state as well as the nature of the transition 
to the paramagnetic state by approximate diagonalization on finite lattices. The recently 
developed density matrix renormalization group method by White J7| and our own exten- 



sion 



26fl to two dimensions allow us to perform calculations on much larger lattices than 
previously possible and with high accuracy 

For L x x L y lattices, when L y > 2L X the difficulty of the calculation depends only weakly 



on L y |26]. By making L y suitably long, spacing between the nearest k y can be made as 
small as we want and thus the close-shell effect can be eliminated. 

Based on diagonalizations of L x x L y lattices with L y = 20, we find the critical hole 
concentrations for the onset of instability in the Nagaoka state to be almost the same for 
L x = 2, 3 and 4. This suggests that the critical hole doping we calculated, S c = 0.22, is close 
to the bulk limit. This value is close to but lower than S r = 0.29 obtained from Edwards trial 



wave function for the case of single spin flip. We show in contrary to previous finding 
2"3fl that a finite region of hole doping exists below 0.22 where the fully ferromagnetic state 
is stable. 

We calculate the ground-state energy of small clusters using the DMRG method, in which 
one reduces degrees of freedom by keeping the eigenstates of the density matrix 0. This 
is in contrast to the conventional real space renormalization group method where the low 
energy eigenstates of the block Hamiltonian are kept. An iterative procedure J7| is used 
to systematically improve the approximation to the density matrix. The DMRG method 
proves to be highly accurate for one dimensional systems. For quantum spin chains, the 
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ground-state energy can be calculated to a high accuracy of 10 6 |7|. When the method 



is applied to the quasi-one dimensional system of several coupled chains pBj, the number 



of states needed to compute the energy to a fixed accuracy grows exponentially with the 



number of chains, but is independent of the length of the chains. It can also be shown ^6 
that the energy calculated in the finite cluster DMRG method always provides a variational 
upper bound to the ground-state energy 

We study the one-band Hubbard model with U = oo on L x x L y square lattices with free 
boundary conditions in both directions. We are restricted to small L x because the accuracy 
of the DMRG method deteriorates at large L x . In this work, the calculations are done on 
strips with L x = 2, 3, 4, 5 and L y = 20. The large value of L y used reduces finite size effects 
due to the k-space shell closing discussed previously. 

Let En(Q,S z ) be the energy calculated for the system with Q holes (the number of 
electrons is iV — Q) and total z-direction spin S z on an L x x L y lattice with N sites. The 
critical hole doping 5 C is determined by comparing E N (Q, S z = 0) with the energy of the 
Nagaoka state, E nag (Q), which is the energy of N — Q spin up electrons on the same lattice. 
(We assume the number of electrons iV — Q is even. When iV — Q is odd, set S z = 1/2.) 
Because of the global SU(2) symmetry of the Hubbard model, the Nagaoka state with total 
spin S = (N — Q)/2 is (2S + l)-fold degenerate. One of these states has S z = 0. Since the 
DMRG method calculates a variational upper bound to the ground-state energy, we have 
En{Q,S z = 0) > E nag (Q), if the ground state is the Nagaoka state. On the other hand, 
if En(Q,S z = 0) < E nag (Q) the ground state is not the Nagaoka state. The smallest hole 
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doping for which this occurs determines the critical doping 5 = Q/N. Since energy computed 
in the DMRG method is a variational upper bound, the critical doping S c we estimated from 
the condition E N (Q, S z = 0) < E nag {Q) is an upper bound to the true S c . 

Since below 5 C the exact E N (Q,S Z = 0) is equal to E nag (Q), the difference between 
the actual E N (Q,S Z = 0) calculated and the corresponding Nagaoka energy provides an 
estimate for the accuracy of our calculations. The accuracy of our calculations in the relevant 
doping region varies from 0.03% for 2 x 20(with M = 52) to 0.5% for 4 x 20 lattices(with 
M = 102)(Fig.[l]), where M is the number of states kept. 

In Fig.0, the energy difference E N (Q, S z = 0) — E nag (Q) between the calculated energy 
and the Nagaoka energy is shown as a function of hole doping (S z = 1/2 if the number 
of electrons is odd). The energies are normalized to E nag (Q). At critical doping 8 C , the 
energy difference changes from positive to negative. The calculated energy En(Q,S z = 0) 
reported here are for the largest number of internal states kept. We have not attempted 
the extrapolation to the infinite- M limit because we are uncertain about the validity of such 
an extrapolation and because at the largest M, En{Q,S z ) gives a nice variational upper 
bound. For 5 < S c , the energy difference is positive and flat. For 5 > 5 C the energy difference 
turns negative abruptly and decreases linearly with 5 — S c (ai least for L x = 2,4). In Fig. 
1(a), the data for 2 x 20 and 2 x 30 are almost indistinguishable from each other indicating 
that L y = 20 is large enough. Since the number of states needed for calculations with fixed 
accuracy grows exponentially with the L x , the errors for 5 x 20 (Fig.|I](d)) are considerably 
larger. Because of this, the energy difference becomes negative at a higher doping. 
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The condition E N (Q,S Z = 0) > E nag (Q) is a necessary condition for the stability of 
the Nagaoka state. It only suggests but does not prove the ground state is ferromagnetic. 
However, for the 2 x 20 and 2 x 30 lattices when hole doping is smaller than 0.22, E^(Q, S z = 
0) — E nag {Q) is as small as 10~ 5 E nag (Q) which strongly suggest that the true E N (Q, S z = 0) 
is in fact equal to E nag (Q). The ferromagnetic state is at least a degenerated ground state. 
The similarity between the data for L x = 3, 4 and L x = 2 suggests that the Nagaoka state is 
stable below about 20 percent doping. Also, the critical hole dopings change very little for 
L x = 2,3,4. This insensitivity indicates that 5 C = 0.22 is close to the bulk limit. 

Near 5 C we also calculated the energy of the lowest state with one spin flipped (S z = 
( N ~Q^ — 1). We expect to achieve higher energy accuracy because the dimensions of the 
Hilbert space is reduced from the S z — case. We have verified that for L x = 2, 3, 4, the 5 C 
inferred from the energy with one spin flipped is the same as the S z = case. 

To investigate the effects of lattice shape anisotropy on the critical doping 5 C , we intro- 
duced hopping anisotropy: t x = 0.5 and t y = 1 on L x x L y lattice with L x = 4 L y = 20. 
Remarkably, the critical doping for this system(Fig.||D is very close to 5 C = 0.22 of the 
isotropic case (when t x — t y — 1 in Fig. 1(c)). This insensitivity to hopping anisotropy gives 
us some confidence that 5 C = 0.22 is close to the bulk limit. 

We now discuss the nature of the ferromagnetic to paramagnetic transition after the 
doping exceeds 5 C . There are two possibilities: (i) the total spin S of the ground state 
changes discontinuously from the maximum N ~® to zero, or (ii) as the hole concentration 
5 exceeds the critical doping, the ground-state total spin reduces gradually to zero as 5 is 
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increased. We can in principle distinguish between these two possibilities by computing 
the ground state energy En(Q,S z ) as a function of S z . In case (i), the energy En(Q,S z ) 
decreases until S z reaches zero. For case (ii), S z stops decreasing at S C (Q) and S C (Q) goes 
to zero when Q is increased. 

Our data is consistent with case (ii) above namely that there exists a doping region 
5 C < 5 < 5 Cl where the ground-state total spin is between S max = — =^ and zero. For 5 > 5 C1 
the ground-state total spin becomes zero. Fig|| shows some representative data. In Fig||(a) 
for 5 = 0.3 on 2 x 30 lattice, the ground-state energy E(Q, S z ) drops quickly with decreasing 
S z until S z = 0.5S max . After that the energy is flat. The total spin of the ground state is 
then S = 0.5S max . The slight increase in the energy from S z = 0.5S max to S z = is due to 
the increased Hilbert space at small S z which makes the calculation less accurate. At 5 = 0.5 
( Fig.|3|(b)), the energy decreases continuously to S z = 0. This implies that the ground state 
has zero total spin. Similar behaviors are observed for L x = 4. We are unable to determine 
8 C1 accurately. But it is close to 0.40. 

We now discuss technical details of the DMRG calculation specific to the infinite-?/ limit. 
A general discussion of DMRG procedures for quasi-one-dimensional systems can be found 



in Ref. |pq| . The chief computational advantage of the infinite-?/ limit over the full Hubbard 
model is the reduced Hilbert dimensions. When expanding a block, we add three states per 
site (empty, spin up and spin down). 



The one dimensional system is used to initialize the environment blocks |7|,p6|. One 
particular difficulty in the infinite-?/ limit is that in one dimension all the spin configurations 
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have exactly the same energy. The total angular momentum is therefore undefined. In the 
quasi-one-dimensional system, however, this degeneracy is lifted and the ground state has 
well defined total angular momentum. One can get around this problem by starting from 
the one dimensional Hubbard model with large U which lifts the degeneracy. 

Typically six iterations are performed for each of several values of M (number of internal 
states kept) starting from small M. To preserve the SU (2) symmetry, we always keep the 
states with the same weight so that the actual number of states kept may be larger than the 
assigned M. 

The programming for the DMRG method is much more complex than Lanczos exact 
diagonalization. Our code for the two dimensional Hubbard model contains over 4000 lines. 
A crucial issue is how to make sure that the computer program is correct. Our computer 
code has passed two non-trivial tests, (i) Hubbard model satisfies a global spin SU(2) 
symmetry. For states with zero total spin, the internal states retained should exhibit the 
SU (2) symmetry. In particular, the states come in 2S + 1 multiplets, i. e. whenever we have 
a state with z-component of spin S z — S we should also find states having identical weight 
(the diagonals of density matrix) with z-direction of spin being S — 1, S — 2, —S + 1, —S. 
(ii) For spin polarized state with S z = ^t^, the computed energy approaches the exact 
answer. 

In conclusion, we have studied the stability of the Nagaoka state in the infinite-f/ Hubbard 
model in two dimensions using the density matrix renormalization group method. We found 
the ferromagnetic state to be stable for a finite doping range near half filling. By computing 



9 



energy upper bounds on L x x 20 lattices with L x up to 5, we have shown that the Nagaoka 
state becomes unstable for hole doping larger than 22 percent. The ground-state total spin 
decreases gradually as the hole doping is increased and becomes zero for more than about 
40 percent doping. 

The work was supported in part by the Office of Naval Research Grant No. N000 14-92- 
J-1340. 
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FIGURES 

FIG. 1. The critical doping S c is obtained by comparing a variational upper bound of energy, 
E, calculated at doping S and S z = (and S z = 1/2 if the number of electrons is odd) with the 
energy of Nagaoka state, E nag . The energy difference, normalized to the Nagaoka energy, is plotted 
as a function of hole doping. The Nagaoka state becomes unstable when the energy difference is 
negative, (a) 2 x 20 and 2 x 30 lattices. The number of states kept in the calculation, M = 52. 
The energy accuracy of the variational bounds are 0.03 percent, (b) 3 x 20 lattice with M=62. (c) 

4 x 20 lattice with M=102. (d) 5 x 20 lattice with M=120. 

FIG. 2. Same as Fig. 1(c) but with anisotropic hoping t x = 0.5 (in the short direction of the 
lattice) and t y = 1. M is 120. 

FIG. 3. The ground-state energy E{S Z ) as a function of S z calculated at M=62. (a) At doping 

5 = 0.3, the energy decreases as the spins are flipped from the Nagaoka state until S z ~ 0.55 max , 
where the S max is the spin of the Nagaoka state. The total spin of the ground state is then close to 
0.5S max - (b) At larger doping 5 = 0.5, the energy decreases continuously. The ground-state total 
spin is zero. 
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